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We compute the entropy of antiferromagnetic quantum spin-s chains in an external magnetic field using 
£Nj ' exact diagonalization and Quantum Monte Carlo simulations. The magnetocaloric effect, i.e., temperature 

variations during adiabatic field changes, can be derived from the isentropes. First, we focus on the example 
of the spin-s = 1 chain and show that one can cool by closing the Haldane gap with a magnetic field. We 
i then move to quantum spin-s chains and demonstrate linear scaling with s close to the saturation field. 

, * ^ | In passing, we propose a new method to compute many low-lying excited states using the Lanczos recursion. 
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£ . 1. Introduction 

O ' 

The magnetocaloric effect, i.e., a temperature change induced by an adiabatic change of an 
external magnetic field was discovered in iron by Warburg in 1881 [1]. This effect also has a long 
history in cooling applications. For example, adiabatic demagnetization of paramagnetic salts was 
the first method to reach temperatures below IK [2]. Among the many experimental but also 
theoretical investigations to follow, we would like to mention in particular the computation of the 
field variation of the entropy in a quantum spin- 1/2 XX Z chain [3]. This investigation is remarkable 
in so far as it is one of the very first numerical computations for an interacting quantum many-body 
system. The relation of the entropy of quantum spin chains in a magnetic field to magnetic cooling 
was also noted [4] and even one exact computation was performed, albeit for a chain with a strong 
Q\ \ easy-axis anisotropy [5]. 

The magnetocaloric effect in quantum spin systems has recently attracted renewed attention. 
The reason for this is two-fold. On the one hand, field-induced quantum phase transitions lead 
to universal responses when the applied field is varied adiabatically [6, 70. On the other hand, it 
was observed that the magnetocaloric effect is enhanced by geometric frustration [9-23] , promising 
improved efficiency in low-temperature cooling applications. In this context, the isotropic spin-1/2 
Heisenberg chain, the isotropic and anisotropic XY chain in a transverse field [10] , and the quantum 
Ising chain in a transverse field [7] were also revisited. One further example of a non-frustrated 
quantum spin chain where the magnetocaloric effect was computed is the ferrimagnetic spin- 1-5/2 
chain [24]. Finally, the adiabatic cooling rate has recently been measured on the spin-1/2 Heisenberg 
chain compound [Cu(/i-C204)(4-aminopyridine)2(H20)]„ [25] and relatedly the magnetocaloric 
effect of the spin-1/2 isotropic Heisenberg and XX Z chains was computed exactly [26]. 

In this paper we will investigate magnetocaloric properties of quantum spin-s Heisenberg chains 
given by the Hamiltonians 

L L 

H = J^2§i-S i+ i-h^2Sf. (1) 

i=l i=l 



1 The presence of precursors should, however, be noted, see, e.g., [4,8]. 
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J > is the antiferromagnetic exchange constant, h an external magnetic held, L the length of the 
chain, and S$ are quantum spin-s operators at site i. We will use periodic boundary conditions, 
i.e., Sl+i — Si. 

In section [5] we will first summarize the methods which we are going to use, namely exact 
diagonalization (ED) and Quantum Monte Carlo (QMC) simulations using the Stochastic Series 
Expansion (SSE) framework. Then we will apply these methods to the s — 1 chain in section [3] 
The spin s = 1 Heisenberg chain is famous for the presence of the so-called 'Haldane' gap at 
h = [27-30]. In particular, we will illustrate with the s = 1 chain that one can exploit the 
presence of a spin gap for cooling by adiabatic magnetization at low magnetic fields. Next, we 
discuss the scaling with spin quantum number s in section [4] In particular, we will show that there 
is a quantum scaling regime at low temperatures close to the saturation field where one has linear 
scaling with s. Finally, in section [5] we will conclude with a summary. 



2. Methods 



2.1. Exact diagonalization 

Thermodynamic quantities can be computed using spectral representations. For example, the 
entropy can be written as [10] 



5 = ^E^ e " £ " /T+lnZ ' 



(2) 



where Z = ~J2 n e~ En ^ T is the partition function and E n are the eigenvalues of the Hamiltonian. 

In order to evaluate spectral representations like ^ , one needs to diagonalize the Hamiltonian. 
If this is done numerically exactly, this is called Exact Diagonalization (ED). First, one should 
perform a symmetry reduction. We have used translations and S^-conservation as well as reflec- 
tions and spin inversion where appropriate. One can then use a library routine to perform a full 
diagonalization. Such an approach is very much in the spirit of classic work [3,31,32]. The main 
technical differences of our computations and [3,31,32] are: (i) we have exploited S'?7(2)-symmetry 
to reconstruct the S z — 1 sector from the S z > 2 sectors and the spin-inversion resolved S z = 
sector, (ii) we have used improved implementations of the diagonalization routines [33,34], (iii) and 
we have substantially more powerful computers at our disposal. The combination of these factors 
enables us to compute full spectra for bigger systems than in the 1960s. 

If we are interested only in low-energy properties, we can use iterative diagonalization algo- 
rithms like the Lanczos method [35-37]. 

The basic Lanczos algorithm [35-37] for a Hcrmitian matrix H proceeds as follows. For a 
given normalized start vector v\ , one defines a sequence of normalized Lanczos vectors Vj via the 
recurrence relations 



0j+l tfj+i 



Hvj-ajVj 



j-i . 



(3) 



0j+i = v ]+ i-Hv 



with the initial conditions (3\ — and vq — 0. The real coefficients oti and fli define the so-called 
Lanczos matrices 

\ 



T„ 



'■• 



Vo 



•. 

'■■ 0n 

n a n J 



(4) 



2 



Magnetocaloric effect in quantum spin-s chains 



The crucial point of the Lanczos algorithm is that the eigenvalues of T n yield a good approximation 
of the extremal eigenvalues of H already for n much smaller than the total Hilbert space dimension. 

The Lanczos algorithm suffers from one practical problem: the recursion relations ([3]) are sup- 
posed to guarantee mutual orthogonality of the vectors Vi, i.e., Vi ■ Vj — for i ^= j. However, 
this fails to be correct when the computations are carried out numerically, leading to undesired 
spurious states, called ghost states [36,37]. Several strategies have been proposed [36,37] in order 
to deal with the ghost problem during the computation of excited eigenvalues. 

Here we suggest another method to compute excited states with controlled accuracy and correct 
multiplicity. We start by performing a fixed number n of Lanczos iterations with a given start vector 
v\. Next we compute the orthogonal transformation Uij, 1 < i, j < n diagonalizing the Lanczos 
matrix ((4|). Then we fix a number m <C n and repeat the Lanczos procedure with the same start 
vector v\. During this second Lanczos pass we construct m vectors 

n 

Hi = Uj t i v j , 1 < i < m . (5) 

3=1 

The second Lanczos pass is needed in order to avoid storing the u> m Lanczos vectors Vj. 

The vectors Ui given by ^ yield approximations to the eigenvectors of H. However, due to the 
ghost problem, there are spurious vectors which need to be eliminated. We perform this in two 
steps. First, using the observation that the vectors u*i are supposed to be mutually orthogonal, we 
can eliminate those Ui which have big projections on the Uj with j < i. Second, we reorthogonalize 
the remaining vectors yielding rh < m orthonormal vectors Ui. Finally, we project the matrix H 
onto this subspace via 

Hij — Ui ■ Hu 3 . (6) 

A full diagonalization of the mxrh matrix Hij and application of the resulting basis transformation 
to the vectors yields rh orthogonal vectors Wi which approximate the eigenvectors of H. The 
eigenvalues Ei associated to the vectors Wi and their accuracy can be estimated according to 

Ei =vS i -Hw i , SEf = w t ■ (H - Eif w, . (7) 

According to our experience, there are neither any missing nor any spurious eigenvalues among 
the converged ones. It should be mentioned that there is no guarantee that the Lanczos procedure 
yields the complete spectrum for a given start vector v\, in particular if one has not performed 
a complete symmetry decomposition. Nevertheless, numerical noise seems to prevent this from 
happening, at least for sufficiently generic start vectors v\ . 

For our procedure one needs to choose the number of Lanczos iterations n and the dimension of 
the subspace m < n in advance and the whole procedure needs to be repeated if it turns out that 
less than the desired number of converged eigenvalues have been obtained. In practice, however, 
one obtains results of a comparable quality for a class of problems where the choice of n and m 
has been adjusted for one representative case. In the examples to be reported below, we have been 
able to obtain about 100 eigenvalues in a given symmetry sector with a relative accuracy of 10 -8 
or better using n = 2000 Lanczos iterations and m = 500 vectors u%. 

2.2. Quantum Monte Carlo 

Since the system sizes accessible by either full diagonalization or Lanczos diagonalization are 
limited, it is desirable to have other methods at our disposal. For non-frustrated spin models like the 
spin-s Heisenberg chain |T]), one can in principle use Quantum Monte Carlo (QMC) simulations. 
There is, however, one problem: we are particularly interested in the entropy which is usually 
obtained with large statistical errors from integrating a Monte-Carlo result for the specific heat. 
In this subsection we will summarize how one can circumvent this problem. 

In order to obtain a QMC estimate of the temperature dependence of the entropy, we employed 
an extended ensemble, broad-histogram method [38-40] within the Stochastic Series Expansion 
(SSE) framework [41-44]. Based on this approach, thermodynamic quantities can be obtained 
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Figure 1. Entropy per site S/L of the s — 1 Heisenberg chain as a function of temperature T 
for h — 3 J. Lines are obtained by ED for L = 14 and 20 sites, symbols show QMC results for 
L = 40 sites. 



over a broad range of temperatures from a single QMC simulation that provides estimates of 
the expansion coefficients g(n) of the system's partition function Z in a high-temperature series 
expansion: 

oc 

Z = Tre-^ = ^ P (n)/3". (8) 

71=0 

Here, (3 = 1/T is the inverse temperature. In terms of the expansion coefficients g(n), the internal 
energy is obtained using 

E = (H) = = -I^n)^- 1 = -hn), (9) 

" n " 

and the free energy from 

F = -±lnZ = -hnJ2g(n)P n . (10) 

n 

Finally, the entropy can be calculated using S = (E — F)/T. 

In the QMC simulation, estimates of the first A coefficients g(n),n = 0, ...,A are obtained by 
performing a random walk in the expansion order n, such as to sample efficiently all expansion 
orders from n — (where g(Q) — (2s + \) L for a spin-s chain of L sites is known exactly) to 
n = A. This is accomplished using an extension of the Wang-Landau flat-histogram sampling 
algorithm [45,46] to the quantum CMSC . £IS detailed previously in Refs. [38-40]. Knowledge of the 
first A coefficients g(n) allows for the calculation of thermodynamic quantities from = down 
to a temperature T m = l//3 m for a system of L sites, where A scales proportional to (3 m x L [38]. 
In particular, for the spin s = 1 Heisenberg chain, we were able to treat systems with up to 
L = 40 sites, accessing temperatures down to T m ps 0.03, which required an already large value of 
A = 5000. The data shown in this contribution (and also those in [14]) were obtained based on the 
original version of the algorithm [38], while recently an improved sampling strategy was proposed, 
based on optimizing the broad histogram Quantum Monte Carlo ensemble [40] . 

2.3. A test case 

We test and compare the aforementioned numerical methods using the example of the s = 1 
Heisenberg chain |T]) with a finite magnetization. Fig. [1] shows the entropy per site S/L as a 



4 



Magnetocaloric effect in quantum spin-s chains 




Figure 2. Isentropes, i.e., adiabatic demagnetization curves of the s = 1 Heisenberg chain as 
obtained from ED and QMC for different chain lengths. The corresponding values of the entropy 
are S/L = 0.05, 0.1, 0.15, 0.2, 0.25, . . . (bottom to top). 

function of temperature T for a fixed magnetic field h = 3 J. For L = 14, we have performed a full 
diagonalization. The corresponding ED curve can therefore be considered as the exact result for 
L = 14. 

A full determination of the spectrum is clearly out of reach already for L — 20. In this case, we 
therefore had to perform a severe truncation to low energies. Thus, the L = 20 ED curve is only 
a low-temperature approximation. In this paper, the precise temperature range is fixed as follows: 
Let £ be the highest energy until which the spectrum is definitely complete. Then we restrict to 
temperatures T < £/5. This choice ensures that missing states are suppressed by a Boltzmann 
factor < exp(— 5) ~ 0.006738. This may seem a small number, but since we are discarding many 
states, extending truncated data to higher temperatures would yield artifacts which are clearly 
visible on the figure(s). 

QMC (denoted by SSE in Fig. Q]) is able to treat the larger system size L = 40. This method, 
however, is best suited for high temperatures. Indeed, at T w 0.05 J one can see deviations caused 
by statistical errors in the SSE result of Fig. [1] which prevented us from going to lower temperatures. 
Finite-size effects can be observed in Fig. [T] in the L — 14 ED curve for T < 0.4 J and in the 
L = 20 ED result. Otherwise it is reassuring that we observe good overall agreement between all 
three methods. 

In passing we note that, for the value of the magnetic field used in Fig.[TJ the low-energy physics 
of the s = 1 chain is described by a Luttingcr liquid (see Ref. [47] and references therein). It is 
well known that the specific heat C of a Luttingcr liquid is linear in T [48]. Due to the relation 
C = T (d S/dT) and because of S(T = 0) = 0, the entropy of a Luttinger liquid is identical to its 
specific heat and in particular also linear in T. Indeed, Fig. [1] is consistent with a linear behavior 
S oc T at sufficiently large L and low T. 

3. Spin s = 1 Heisenberg chain 

In this section we will discuss the entropy of the s = 1 Heisenberg chain in more detail. 

First, we show results for the value of the entropy per site S/L as a function of h and T in 
Fig.[5J The ED curves have been obtained by computing the entropy S on a mesh in the /i-T-planc 
and determining the constant entropy curves, i.e., the isentropes from this data. This leads to some 
discretization artifacts at the cusps in the the S/L = 0.05 curves (the ED curves at the lowest 
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Figure 3. Entropy per spin of the s = 1 Heisenberg chain at zero magnetic field h — as a 
function of temperature. Note that the temperature axis is scaled as 1/T and the entropy axis 
logarithmically. The dashed line is a fit to flT2} with A = 0.577 and A = 0.4105 J. 

temperature) in Fig. [2] Conversely, the QMC had problems to resolve the large low-temperature 
entropy around the saturation field h aa ,t. = 4 J. Accordingly, the L — 40 SSE data points with 
S/L = 0.05 are missing at h/ J = 3.7 and 4. From our QMC simulations we can just conclude that 
for an L = 40 chain these two points are in the region T < 0.067 J and T < 0.085 J for h = 3.7 J 
and 4 J, respectively. 

The ED curves with S/L = 0.05 in Fig. fj] have marked finite-size wiggles. Also at higher 
temperatures (entropies), finite-size effects can be observed. There are, however, two regions where 
finite-size effects are evidently small. One is the low-temperature regime for fields h smaller than 
the Haldane gap [29, 30] 



Also in the gapped high-field regime above the saturation field h>AJ one observes essentially no 
finite-size effects at any temperature. This renders ED the method of choice at high fields h > h sat ., 
in particular at low temperatures, while QMC is preferable otherwise because bigger system sizes 
L are accessible. 

One can read off from Fig. [2] pronounced temperature changes in two regimes. Firstly, one can 
obtain cooling by adiabatic demagnetization from a high magnetic field h > h aat . as one lowers 
the field to the saturation field. Secondly, at low temperatures, one can also cool by adiabatic 
magnetization from h = to h = A. This demonstrates that one can use the gap-closing at a 
generic field-induced quantum phase transition for cooling purposes. In the regime 0.41 05 J ~ A < 
h < hsat, = 4 J, the spectrum is gapless (compare Fig. Q] and the related discussion in section l2~3l) . 
Accordingly, one observes only small temperature changes induced by adiabatic (de)magnetization 
in this field range. 

We now focus on the low-field regime. Firstly, we consider the behavior of the entropy at h = 
which is shown in Fig. [31 One observes that finite-size effects are small for L > 14: the L = 14 
curve is almost within the error margins of the L — 40 QMC curve even at low temperatures. This 
can be attributed to a finite correlation length £ which is evidently sufficiently smaller than 14 for 
all temperatures at h = 0. Indeed, the correlation length of the s = 1 Heisenberg chain is known 
to be £ w 6 at h = and T = [29,30]. 

Because of the presence of a gap A at h = 0, we expect the entropy to be exponentially 
activated as a function of temperature. Accordingly, the low-temperature asymptotic behavior at 



A « 0.4105 J. 



(11) 
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Figure 4. Entropy per spin of the s = 1 Heisenberg chain at h = 0.4105 J, corresponding to the 
Haldane gap. Note the doubly logarithmic scale. The dashed line shows the square-root behavior 
(fj"3l> with B = 0.1323. 



h = should follow the form 



«4e' 



-A/T 



(12) 



/i=0 



If we fix the gap to the known value A = 0.4105 J [29,30], the only free parameter in the low- 
temperature asymptotic behavior (fT2]) is the prefactor A. A fit of the L — 40 QMC results yields 
A = 0.577(1). The dashed line in Fig. [3]shows that the formula (fT2|) describes the low-temperature 
regime quite well with the given parameters. 

Secondly, we consider a magnetic field exactly equal to the Haldane gap h = 0.4105 J. Fig. [4] 
shows a doubly-logarithmic plot of the entropy S as a function of temperature T for this value of 
the magnetic field. Since we are now sitting at a field-induced quantum phase transition, we expect 
a large number of low- lying states and relatedly an infinite correlation length £ = oo. On a finite 
system, these low-lying states will be pushed to higher energies. Indeed, we observe pronounced 
finite-size effects at low temperatures in Fig. |4] 

At a quantum phase transition in one dimension which preserves a U (l)-symmetry, the entropy 
S is expected [4,6,7,10] to vary asymptotically as a square root of temperature: 



h=A 



= 8,/ 7 . 



(13) 



As before, the amplitude B is the only free parameter. A fit to the L = 40 QMC data with 
T < 0.07 J yields B = 0.1323(6). The dashed line in Fig. H shows that the L — AO QMC data is 
indeed consistent with the asymptotic square root (flU)) in the range 0.03 <T/J< 0.07. However, 
because of finite-size effects, in this case we really need to go to L = 40 to recover the asymptotic 
behavior. Furthermore, we need to restrict to T < 0.07 J since one can see from Fig. [4]that a pure 
power law is no longer a good description for T > 0.07 J. 

Comparison of (|12p and (|13[) shows that one can reach exponentially small temperatures Tf 
at h = A by adiabatic magnetization of a spin-1 Heisenberg chain from an initial temperature T; 
at h = 0. For an initial temperature T\ < 0.1 J, the final temperature Tf can be quantitatively 
estimated from the following combination of (fT2|) and (| 13|) 



Tj__A*_ 
J B 2 

with the parameters A, A, and B given above. 



-2 A /Ti 



(14) 
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Figure 5. Isentropes of L = 20 Heisenberg chains as obtained by ED for different s. Note that 
a truncation of the energy spectrum has been used for s > 1/2. The values of the entropy per 
site are S/L = 0.05, 0.01, 0.15, 0.2, and 0.25 (bottom to top). 

4. Scaling with spin quantum number s 

We will now turn to the high-field region and consider in particular scaling with the spin 
quantum number s close to the saturation field. 

It is straightforward to compute the excitation energy e of a single flipped spin propagating 
above a ferromagnetically polarized background (see, e.g., [49]): 

s{k) = 2 J s (cosfc- 1) +h. (15) 

The minimum of this one-magnon dispersion is located at k = tt. The saturation field h sat . of 
the spin-s Heisenberg chain ([T]) is given by the one-magnon instability [50]. Accordingly, it can be 
determined by inserting the condition e{n) = into (| 1 5[) : 

h sat . = 4Js. (16) 

The fact that the single-particle energy (fT5")) scales linearly with s suggests to scale all energies 
close to the saturation field and at low temperatures linearly with s. In order to test this scenario, 
we have first computed the isentropes for a fixed chain length L = 20 and s < 2. The s = 1/2 
curves are obtained by a full diagonalization of the Hamiltonian [10] while we have performed 
additional ED computations for s = 1, 3/2, and 2 using the truncation procedure described in 
section [2] Fig. [5] shows the resulting isentropes for S/L = 0.05, 0.01, 0.15, 0.2, and 0.25 close to 
the saturation field (fTB)) . We observe that at low temperatures, scaling both h and T by s leads to 
a nice collapse of the isentropes around and in particular above the saturation field. 

Further details can be read off from the temperature scans at a fixed magnetic field shown in 
Fig. [6] In this figure, we have chosen to show the specific heat C as a representative for energy- 
related quantified and the magnetic susceptibility \ as a representative for magnetic quantities. 
The top row of Fig. [5] corresponds to a magnetic field h = 1.1 h sa t.- Here we observe a nice collapse 
with s for T / (J s) < 0.15. The middle row of Fig. [5] corresponds to a magnetic field exactly equal to 
the saturation field. Here we observe a good scaling collapse for T / (J s) < 0.04. Finally, the bottom 
row of Fig. [6] corresponds to a magnetic field h = 0.9/i sa t.. In this case, one expects the scaling 
region to be pushed to even lower temperatures. However, for h < h sgut . the finite size L — 20 leads 

2 The specific heat C = T (d S/dT) is equivalent to S since in the present case S(0)/L — » for L — » oo. 
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Figure 6. Specific heat per site C/L (left column) and magnetic susceptibility per site x/L 
(right column) as obtained by ED for different s. Note that a truncation of the energy spectrum 
has been used for s > 1/2. The values of the magnetic field are chosen somewhat above the 
saturation field (top row), exactly at the the saturation field (middle row), and somewhat below 
the saturation field (bottom row). 



to artifacts in the low-temperature behavior. Therefore it is difficult to make definite statements 
for this case. 

To summarize this section, we have provided evidence of linear scaling with s for thermodynamic 
quantities of spin-s Heisenberg chains close to the saturation field. Note that this is very different 
from the quadratic scaling, e.g. of T with s 2 [23] or with s (s + 1) [51,52], which is needed to 
approach the classical limit. 
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5. Conclusions and outlook 

In this paper we have illustrated the numerical computation of thermodynamic quantities and 
in particular the entropy for spin-s Heisenberg chains. From a technical point of view, we have 
described in section 12.11 how to perform a reliable computation of a large number of low-lying 
states using the Lanczos method and in section 12.21 how to compute the entropy directly by a 
QMC simulation. 

In section [3] we have then focused on the s = 1 Heisenberg chain and shown in particular that 
one can cool with an adiabatic magnetization process during which the Haldane gap A is closed. 
Many previous investigations (e.g., [9, 10, 14, 18, 19, 23]) have focused on the saturation field. The 
reason for this is that the saturation field is a field-induced quantum phase transition at a known 
value of the magnetic field which also gives rise to technical simplifications. However, the scenario 
of quantum phase transitions [6, 7] is universal and not restricted to the saturation field. Indeed, 
Fig. 2] is consistent with the same universal square-root behavior of the entropy S at h = A in 
the spin-1 Heisenberg chain as observed previously for s — 1/2 chains exactly at the saturation 
field [4, 10]. This may also be important from an experimental point of view since a possible spin 
gap may be accessible by laboratory magnetic fields even if the saturation field is out of reach. In 
fact, cooling by adiabatic magnetization when closing a spin gap has presumably been indirectly 
observed in pulse-field magnetization experiments on SrCu2 (603)2 [53]. 

In section [4] we then moved to the saturation field and investigated scaling with the spin 
quantum number s. In contradistinction to the classical scaling regime where temperature should 
scale quadratically with s [23,51,52], there is a quantum scaling regime around the saturation 
field where one observes a collapse using the linearly scaled parameters h/s and T/s. Also in 
higher dimensions the one-magnon dispersion typically scales with s. By the same arguments as 
in section [4j we therefore expect linear scaling with s at any continuous transition to saturation. 
This expectation could be tested numerically with the methods of the present paper. 
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